######################################
# DEMOCRATIC PEACE SURVEY CODE

# clear work space
rm(list = ls())

# print date and time of run
print("Starting date and time")
Sys.time()

# load packages and functions
R_folder <- "INSERT HERE"
source(paste0(R_folder, "survey_functions.R"))

# display the output
output_print <- TRUE
# return the main results
return_main <- TRUE

# images directory
images_directory <- "INSERT HERE"

# grayscale for the graphics or not
#print_colormodel <- "srgb" # in color
print_colormodel <- "grey"

# load data
data_folder <- "INSERT HERE"
load(paste0(data_folder, "democratic_peace.RData"))

##########################################3
# Treatment Assignments Variables
# Democracy Vignette
d$Z <- ifelse(d$d03==1 | d$d08==1 |
                d$d09==1 | d$d10==1 | d$d11==1 | d$d13==1 | d$d14==1, 1, 0)
d$Z[d$d02==1 | d$d04==1 | d$d05==1 | d$d06==1 |
      d$d07==1 | d$d12==1 | d$d15==1] <- 0 
d$Z <- factor(x = d$Z, labels = c("Non-democracy", "Democracy"))
# Vignette Types
# Basic 
d$V[d$d02==1 | d$d03==1] <- 1
# Covariate Control
d$V[d$d04==1 | d$d05==1 | d$d06==1 |
      d$d07==1 | d$d08==1 | d$d09==1 |
      d$d10==1 | d$d11==1] <- 2
# ENE
d$V[d$d12==1 | d$d13==1 | d$d14==1 | d$d15==1] <- 3
d$ENE <- ifelse(d$d12==1 | d$d13==1, "Fragile Democracy", NA)
d$ENE[d$d14==1 | d$d15==1] <- "Fragile Dictatorship"
d$ENE <- factor(d$ENE)
# make a factor for d$V
d$V <- factor(x = d$V, labels = c("Basic", "Covariate Control", "ENE"))

# remove the ones missing treatment assignment 
d <- d[!is.na(d$V),]

# size of data set
if (output_print) {
  dim(d)
}

# check whether treatment assignment affects completing the survey in full
d$complete <- !d$workerid_missing
robustse(lm(complete ~ Z, data = d))
table(d$complete, d$V)
fisher.test(table(d$complete, d$V))

############################################
# Placebo Test: Most Likely Countries
d$ocountry <- ifelse(d$ocou_1 == "", d$ocou_2, d$ocou_1) %>% tolower()
country_list <- paste(d$ocountry, collapse = " ")
country_list <- strsplit(x = country_list, split = " ")[[1]]
table_country <- as.data.frame(table(country_list))
table_country <- table_country[order(table_country$Freq),]
# North Korea
d$north_korea <- grepl(pattern = "north korea", d$ocountry)
# Iran
d$iran <- grepl(pattern = "iran", d$ocountry)
# Iraq
d$iraq <- grepl(pattern = "iraq", d$ocountry)
# Saudi
d$saudi_arabia <- grepl(pattern = "saudi arabia|saudi|arabia", d$ocountry)
# Russia
d$russia <- grepl(pattern = "russia", d$ocountry)
# China
d$china <- grepl(pattern = "china", d$ocountry)
# Pakistan 
d$pakistan <- grepl(pattern = "pakistan", d$ocountry)
# Syria
d$syria <- grepl(pattern = "syria", d$ocountry)
# South Korea
d$south_korea <- grepl(pattern = "south korea", d$ocountry)
# Egypt
d$egypt <- grepl(pattern = "egypt", d$ocountry)

c_summary <- d %>% group_by(V, Z) %>%
  dplyr::summarise(north_korea = sum(north_korea, na.rm = TRUE),
            iran = sum(iran, na.rm = TRUE),
            iraq = sum(iraq, na.rm = TRUE),
            saudi_arabia = sum(saudi_arabia, na.rm = TRUE),
            russia = sum(russia, na.rm = TRUE),
            china = sum(china, na.rm = TRUE),
            pakistan = sum(pakistan, na.rm = TRUE),
            syria = sum(syria, na.rm = TRUE),
            south_korea = sum(south_korea, na.rm = TRUE),
            egypt = sum(egypt, na.rm = TRUE))
c_summary <- reshape2::melt(as.data.frame(c_summary))
country_list <- c("North Korea", "Iran", "Iraq", "Saudi Arabia", "Russia",
                  "China", "Pakistan", "Syria", "South Korea", "Egypt")
c_summary$variable <- relab(old.var = c_summary$variable, 
                            old.labs = unique(c_summary$variable), 
                            new.labs = country_list, text = TRUE)
c_summary$variable <- factor(c_summary$variable, 
                             levels = c("South Korea", "Egypt",
                                        "Saudi Arabia", "Syria",
                                        "Pakistan", "Russia", "China",
                                        "Iraq", "North Korea", "Iran"))
c_summary$Z <- factor(as.character(c_summary$Z), rev(levels(c_summary$Z)))
ggplot(data = c_summary, mapping = aes(x = variable, y = value, fill = Z)) + 
  geom_bar(stat = "identity", position = "dodge") +
  facet_wrap(~ V, ncol = 3) + theme_bb() + coord_flip() +
  xlab("10 Most Frequently Mentioned Countries") +  
  ylab("Number of Respondents Who Mentioned the Country") +
  scale_fill_manual(values = c("steelblue4", "orange1"),
                    name = "Treatment Assignment Regime Type",
                    guide = guide_legend(reverse=TRUE))
ggsave(paste0(images_directory, "Figure_22.pdf"), 
       width = 8, height = 6, dpi = 300, colormodel = print_colormodel)

############################################h
# Placebo Test: Regions of the World

# 1) Dimension Reduction
d$NoAm <- region.lab(region1 = d$regions_1_1_Group, region2 = d$regions_2_1_Group)
d$WeEu <- region.lab(d$regions_1_4_Group, d$regions_2_4_Group)
d$MiEa <- region.lab(d$regions_1_6_Group, d$regions_2_6_Group)
d$SuAf <- region.lab(d$regions_1_7_Group, d$regions_2_7_Group)
d$EaAs <- region.lab(d$regions_1_8_Group, d$regions_2_8_Group)
d$CeAs <- region.lab(d$regions_1_9_Group, d$regions_2_9_Group)
d$region <- 1*d$NoAm + 1*d$WeEu - 1*d$MiEa - 1*d$SuAf
summary(d$region)
# graphics
region_summary <- d %>% group_by(V, Z) %>%
  dplyr::summarise(s_NoAm = mean(NoAm, na.rm = TRUE),
            s_WeEu = mean(WeEu, na.rm = TRUE),
            s_MiEa = mean(MiEa, na.rm = TRUE),
            s_SuAf = mean(SuAf, na.rm = TRUE),
            s_EaAs = mean(EaAs, na.rm = TRUE),
            s_CeAs = mean(CeAs, na.rm = TRUE))

region_summary <- reshape2::melt(as.data.frame(region_summary))
region_summary$variable <- relab(old.var = region_summary$variable, 
                                 old.labs = levels(region_summary$variable), 
                                 new.labs = c("North America", "Western Europe",
                                              "Middle East/North Africa", "Sub-Saharan Africa",
                                              "East Asia", "Central Asia"), TRUE)
region_labs <- c("North America", "Western Europe",
                 "Central Asia", "East Asia", 
                 "Sub-Saharan Africa", "Middle East/North Africa")
region_summary$variable <- factor(region_summary$variable,
                                  levels = region_labs)
region_summary$Z <- factor(as.character(region_summary$Z),
                           rev(levels(region_summary$Z)))
ggplot(region_summary, mapping = aes(y = value, x = variable, fill = Z)) + 
  geom_bar(stat = "identity", position = "dodge") + 
  facet_wrap( ~ V, ncol=2) + coord_flip() + theme_bb() + 
  ylab("Proportion of Respondents Who Selected It as a Most Likely Region") + 
  xlab("Geographic Regions") +
  scale_fill_manual(values = c("steelblue4", "orange1"), name = "Treatment Assignment Regime Type") + ylim(0, 1) +
  guides(fill = guide_legend(reverse = TRUE))
ggsave(paste0(images_directory, "Figure_12.pdf"), 
       width = 8, height = 5.5, dpi = 300, colormodel = print_colormodel)
# standardized
d$region.s <- vscale(var = d$region, vig = d$V)
# regression: standardized
region.s <- myreg(Y = d$region.s, Z = d$Z, V = d$V)
colnames(region.s) <- c("coef", "SE")
# regression: non-standardized
region.n <- myreg(Y = d$region, Z = d$Z, V = d$V)
colnames(region.n) <- c("coef", "SE")
# regression with dummy variables and interactions
region.1 <- robustse(lm(formula = region ~ Z + V + Z:V, data = d))
region.2 <- robustse(lm(formula = region.s ~ Z + V + Z:V, data = d))
if (output_print) {
  print("Placebo Test: Regions of the World")
  print(region.s)
  print(region.n)
  print(region.1)
  print(region.2)
}

############################################
# Placebo Test: GDP per Capita
d$gdp <- ifelse(is.na(d$gdp_1), d$gdp_2, d$gdp_1)
d$gdp.o <- d$gdp-1
summary(d$gdp)
# relabel 
# real-world median values for each interval
qog <- read.dta(paste0(data_folder, "qog_bas_ts_jan15.dta"))
r_gdp <- qog$gle_cgdpc[qog$year==2005]
d$gdp <- relab(old.var = d$gdp, old.labs = 1:7, 
               new.labs = c(median(r_gdp[r_gdp < 500], TRUE), 
                            median(r_gdp[r_gdp > 500 & r_gdp < 1001], TRUE),
                            median(r_gdp[r_gdp > 1000 & r_gdp < 5001], TRUE),
                            median(r_gdp[r_gdp > 5000 & r_gdp < 10001], TRUE),
                            median(r_gdp[r_gdp > 10000 & r_gdp < 20001], TRUE),
                            median(r_gdp[r_gdp > 20000 & r_gdp < 40001], TRUE),
                            median(r_gdp[r_gdp > 40000], TRUE)))
d$gdp <- log(d$gdp)
# standardize
d$gdp.s <- vscale(var = d$gdp, vig = d$V)
d$gdp.s.o <- vscale(var = d$gdp.o, vig = d$V)
# regression: standardized log dollar
gdp.s <- myreg(Y = d$gdp.s, Z = d$Z, V = d$V)
colnames(gdp.s) <- c("coef", "SE")
# regression: standardized ordinal scale
gdp.s.o <- myreg(Y = d$gdp.s.o, Z = d$Z, V = d$V)
colnames(gdp.s.o) <- c("coef", "SE")
# regression: non-standardized
gdp.n <- myreg(Y = d$gdp, Z = d$Z, V = d$V)
colnames(gdp.n) <- c("coef", "SE")
# regression with dummy variables and interactions
gdp.1 <- robustse(lm(formula = gdp ~ Z + V + Z:V, data = d))
gdp.2 <- robustse(lm(formula = gdp.s ~ Z + V + Z:V, data = d))
if (output_print) {
  print("Placebo Test: GDP per Capita")
  print(gdp.s)
  print(gdp.s.o)
  print(gdp.n)
  print(gdp.1)
  print(gdp.2)
}

############################################
# Placebo Test: Religion
d$religion <- ifelse(is.na(d$religion_1), d$religion_2, d$religion_1)
summary(d$religion)
# relabel 
probint <- c(mean(0:20), mean(21:40), mean(41:60), mean(61:80), mean(81:100))
likelylab <- c("Very Unlikely", "Unlikely", "Odds About Even", 
               "Likely", "Very Likely")
d$religion <- relab(old.var = d$religion, old.labs = 1:5, 
                    new.labs = probint)
# standardize
d$religion.s <- vscale(var = d$religion, vig = d$V)
# regression: standardized 
religion.s <- myreg(Y = d$religion.s, Z = d$Z, V = d$V)
colnames(religion.s) <- c("coef", "SE")
# regression: non-standardized
religion.n <- myreg(Y = d$religion, Z = d$Z, V = d$V)
colnames(religion.n) <- c("coef", "SE")
# regression with dummy variables and interactions
religion.1 <- robustse(lm(formula = religion ~ Z + V + Z:V, data = d))
religion.2 <- robustse(lm(formula = religion.s ~ Z + V + Z:V, data = d))
if (output_print) {
  print("Placebo Test: Religion")
  print(religion.s)
  print(religion.n)
  print(religion.1)
  print(religion.2)
}

############################################
# Placebo Test: Oil Reserves
d$oil <- ifelse(is.na(d$oil_1), d$oil_2, d$oil_1)
summary(d$oil)
# relabel 
d$oil <- relab(old.var = d$oil, old.labs = 1:5, 
               new.labs = rev(probint))
d$oil.o <- relab(old.var = ifelse(is.na(d$oil_1), d$oil_2, d$oil_1), 
                 old.labs = 1:5, 
                 new.labs = probint)
# standardize
d$oil.s <- vscale(var = d$oil, vig = d$V)
# regression: standardized 
oil.s <- myreg(Y = d$oil.s, Z = d$Z, V = d$V)
colnames(oil.s) <- c("coef", "SE")
# regression: non-standardized
oil.n <- myreg(Y = d$oil, Z = d$Z, V = d$V)
colnames(oil.n) <- c("coef", "SE")
# regression with dummy variables and interactions
oil.1 <- robustse(lm(formula = oil ~ Z + V + Z:V, data = d))
oil.2 <- robustse(lm(formula = oil.s ~ Z + V + Z:V, data = d))
if (output_print) {
  print("Placebo Test: Oil Reserves")
  print(oil.s)
  print(oil.n)
  print(oil.1)
  print(oil.2)
}

############################################
# Placebo Test: White
d$white <- ifelse(is.na(d$white_1), d$white_2, d$white_1)
summary(d$white)
# relabel 
probint <- c(mean(0:20), mean(21:40), mean(41:60), mean(61:80), mean(81:100))
likelylab <- c("Very Unlikely", "Unlikely", "Odds About Even", 
               "Likely", "Very Likely")
d$white <- relab(old.var = d$white, old.labs = 1:5, 
                 new.labs = probint)
# standardize
d$white.s <- vscale(var = d$white, vig = d$V)
# regression: standardized 
white.s <- myreg(Y = d$white.s, Z = d$Z, V = d$V)
colnames(white.s) <- c("coef", "SE")
# regression: non-standardized
white.n <- myreg(Y = d$white, Z = d$Z, V = d$V)
colnames(white.n) <- c("coef", "SE")
# regression with dummy variables and interactions
white.1 <- robustse(lm(formula = white ~ Z + V + Z:V, data = d))
white.2 <- robustse(lm(formula = white.s ~ Z + V + Z:V, data = d))
if (output_print) {
  print("Placebo Test: White")
  print(white.s)
  print(white.n)
  print(white.1)
  print(white.2)
}

############################################
# Placebo Test: Military Spending
d$force <- ifelse(is.na(d$force_1), d$force_2, d$force_1)
# ordinal coding
d$force.o <- d$force
d$force.o <- relab(old.var = d$force.o, old.labs = 1:5, new.labs = 0:4)
# relabel 
# real-world median values in each interval
nmc <- read.csv(paste0(data_folder, "NMC_v4_0.csv"))
r_exp <- nmc$milex[nmc$year == 2005]/1000
d$force <- relab(old.var = d$force, old.labs = 1:5, 
                 new.labs = c(median(r_exp[r_exp >= 0 & r_exp <= 30], TRUE),
                              median(r_exp[r_exp > 30 & r_exp <= 120], TRUE),
                              median(r_exp[r_exp > 120 & r_exp <= 600], TRUE),
                              median(r_exp[r_exp > 600 & r_exp <= 3500], TRUE),
                              median(r_exp[r_exp > 3500], TRUE)))
d$force <- log(d$force)
# standardize
d$force.s <- vscale(var = d$force, vig = d$V)
d$force.s.o <- vscale(var = d$force.o, vig = d$V)
# regression: standardized log million dollar scale
force.s <- myreg(Y = d$force.s, Z = d$Z, V = d$V)
colnames(force.s) <- c("coef", "SE")
# regression: standardized ordinal scale
force.s.o <- myreg(Y = d$force.s.o, Z = d$Z, V = d$V)
colnames(force.s.o) <- c("coef", "SE")
# regression: non-standardized
force.n <- myreg(Y = d$force, Z = d$Z, V = d$V)
colnames(force.n) <- c("coef", "SE")
# regression: ordinal coding
force.o <- myreg(Y = d$force.o, Z = d$Z, V = d$V)
colnames(force.o) <- c("coef", "SE")
# regression with dummy variables and interactions
force.1 <- robustse(lm(formula = force ~ Z + V + Z:V, data = d))
force.2 <- robustse(lm(formula = force.s ~ Z + V + Z:V, data = d))
if (output_print) {
  print("Placebo Test: Military Capability")
  print(force.s)
  print(force.n)
  print(force.o)
  print(force.1)
  print(force.2)
}

############################################
# Placebo Test: Military Alliance
d$allies <- ifelse(is.na(d$alliance_1), d$alliance_2, d$alliance_1)
summary(d$allies)
# relabel 
d$allies <- relab(old.var = d$allies, old.labs = 1:5, 
                  new.labs = probint)
# standardize
d$allies.s <- vscale(var = d$allies, vig = d$V)
# regression: standardized 
allies.s <- myreg(Y = d$allies.s, Z = d$Z, V = d$V)
colnames(allies.s) <- c("coef", "SE")
# regression: non-standardized
allies.n <- myreg(Y = d$allies, Z = d$Z, V = d$V)
colnames(allies.n) <- c("coef", "SE")
# regression with dummy variables and interactions
allies.1 <- robustse(lm(formula = allies ~ Z + V + Z:V, data = d))
allies.2 <- robustse(lm(formula = allies.s ~ Z + V + Z:V, data = d))
if (output_print) {
  print("Placebo Test: Military Alliance")
  print(allies.s)
  print(allies.n)
  print(allies.1)
  print(allies.2)
}

############################################
# Placebo Test: Trade with the U.S.
d$trade <- ifelse(is.na(d$trade_1), d$trade_2, d$trade_1)
summary(d$trade)
# ordinal coding
d$trade.o <- d$trade
d$trade.o <- relab(old.var = d$trade.o, old.labs = 1:5, new.labs = 0:4)
# relabel 
# get real-world median values for each interval
cow_trade <- read.csv(paste0(data_folder, "dyadic_trade_3.0.csv"))
cow_trade <- cow_trade[cow_trade$importer1=="United States of America" & 
                         cow_trade$year==2005,]
r_trade <- cow_trade$flow1 + cow_trade$flow2
d$trade <- relab(old.var = d$trade, old.labs = 1:5, 
                 new.labs = c(median(r_trade[r_trade >= 0 & r_trade <= 100], TRUE),
                              median(r_trade[r_trade > 100 & r_trade <= 350], TRUE),
                              median(r_trade[r_trade > 350 & r_trade <= 1500], TRUE),
                              median(r_trade[r_trade > 1500 & r_trade <= 10000], TRUE),
                              median(r_trade[r_trade > 10000], TRUE)))
d$trade <- log(d$trade)
# standardize
d$trade.s <- vscale(var = d$trade, vig = d$V)
d$trade.s.o <- vscale(var = d$trade.o, vig = d$V)
# regression: standardized log million dollars
trade.s <- myreg(Y = d$trade.s, Z = d$Z, V = d$V)
colnames(trade.s) <- c("coef", "SE")
# regression: standardized ordinal scale
trade.s.o <- myreg(Y = d$trade.s.o, Z = d$Z, V = d$V)
colnames(trade.s.o) <- c("coef", "SE")
# regression: non-standardized
trade.n <- myreg(Y = d$trade, Z = d$Z, V = d$V)
colnames(trade.n) <- c("coef", "SE")
# regression: ordinal coding
trade.o <- myreg(Y = d$trade.o, Z = d$Z, V = d$V)
colnames(trade.o) <- c("coef", "SE")
# regression with dummy variables and interactions
trade.1 <- robustse(lm(formula = trade ~ Z + V + Z:V, data = d))
trade.2 <- robustse(lm(formula = trade.s ~ Z + V + Z:V, data = d))
if (output_print) {
  print("Placebo Test: Trade with the U.S.")
  print(trade.s)
  print(trade.n)
  print(trade.o)
  print(trade.1)
  print(trade.2)
}

############################################
# Placebo Test: Joint Military Exercise
d$exercise <- ifelse(is.na(d$exercise_1), d$exercise_2, d$exercise_1)
summary(d$exercise)
# relabel 
d$exercise <- relab(old.var = d$exercise, old.labs = 1:5, 
                    new.labs = probint)
# standardize
d$exercise.s <- vscale(var = d$exercise, vig = d$V)
# regression: standardized 
exercise.s <- myreg(Y = d$exercise.s, Z = d$Z, V = d$V)
colnames(exercise.s)  <- c("coef", "SE")
# regression: non-standardized
exercise.n <- myreg(Y = d$exercise, Z = d$Z, V = d$V)
colnames(exercise.n) <- c("coef", "SE")
# regression with dummy variables and interactions
exercise.1 <- robustse(lm(formula = exercise ~ Z + V + Z:V, data = d))
exercise.2 <- robustse(lm(formula = exercise.s ~ Z + V + Z:V, data = d))
if (output_print) {
  print("Placebo Test: Joint Military Exercise")
  print(exercise.s)
  print(exercise.n)
  print(exercise.1)
  print(exercise.2)
}

############################################
# Placebo Test: Foreign Direct Investment
sizelab <- c("Very Low", "Low", "Medium", "High", "Very High")
d$invest <- ifelse(is.na(d$invest_1), d$invest_2, d$invest_1)
summary(d$invest)
# relabel 
d$invest <- relab(old.var = d$invest, old.labs = 1:5, 
                  new.labs = 0:4)
# standardize
d$invest.s <- vscale(var = d$invest, vig = d$V)
# regression: standardized 
invest.s <- myreg(Y = d$invest.s, Z = d$Z, V = d$V)
colnames(invest.s) <- c("coef", "SE")
# regression: non-standardized
invest.n <- myreg(Y = d$invest, Z = d$Z, V = d$V)
colnames(invest.n) <- c("coef", "SE")
# regression with dummy variables and interactions
invest.1 <- robustse(lm(formula = invest ~ Z + V + Z:V, data = d))
invest.2 <- robustse(lm(formula = invest.s ~ Z + V + Z:V, data = d))
if (output_print) {
  print("Placebo Test: Foreign Direct Investment")
  print(invest.s)
  print(invest.n)
  print(invest.1)
  print(invest.2)
}

#######################################
# Make coefficient Plot
# standardized: log dollar scales for GDP per capita, 
# military spending, and trade
res.s <- data.frame(rbind(region.s, gdp.s, religion.s, white.s,
                          oil.s, allies.s, trade.s, 
                          exercise.s, invest.s, force.s))
names(res.s) <- c("coef", "se")
res.s$v <- rep(levels(d$V), 10)
# vignette labels
res.s$v2 <- factor(res.s$v, levels=levels(factor(res.s$v)))
res.s$v3 <- factor(res.s$v, levels=rev(levels(factor(res.s$v))))
p.labs <- c("Most Likely Region", 
            "GDP per Capita",
            "Likelihood of Being Majority Christian",
            "Likelihood of Being Majority White",
            "Likelihood of Not Having Large Oil Reserves",
            "Likelihood of Military Alliance with U.S.*",
            "Trade with U.S.*",
            "Likelihood of Joint Military Exercise with U.S.**",
            "Level of Investment in U.S. Businesses**",
            "Military Spending*")
p.labs.func <- function(x, times=3) {rep(p.labs[x],times)}
res.s$placebo <- unlist(lapply(1:length(p.labs), p.labs.func))
res.s$placebo <- factor(res.s$placebo, levels = p.labs)

p.labs_l <- c("A: Most Likely Region", 
              "B: GDP per Capita",
              "C: Likelihood of Being Majority Christian",
              "D: Likelihood of Being Majority White",
              "E: Likelihood of Not Having Large Oil Reserves",
              "F: Likelihood of Military Alliance with U.S.*",
              "G: Trade with U.S.*",
              "H: Likelihood of Joint Military Exercise with U.S.**",
              "I: Level of Investment in U.S. Businesses**",
              "J: Military Spending*")

# standardized: ordinal scales for GDP per capita, 
# military spending, and trade 
res.s.o <- res.s
res.s.o[,1:2] <- data.frame(rbind(region.s, gdp.s.o, religion.s, white.s,
                                  oil.s, allies.s, trade.s.o, 
                                  exercise.s, invest.s, force.s.o))
coef_top <- ggplot(res.s.o[1:15,], aes(x=coef, y=v3, shape=v2, color=v2, alpha = v2))

res.s$placebo <- rep(p.labs_l, each = 3)
res.s.o$placebo <- rep(p.labs_l, each = 3)

# coefficient plot with military spending (standardized)
coef_plot_function(start_ggplot = ggplot(res.s, aes(x=coef, y=v3, 
                                                    shape=v2, color=v2, alpha = v2)), 
                   file_name = paste0("Figure_02"),
                   main_title_name = NULL,
                   x_title_name = "Standardized Difference (Dem - NonDem)",
                   dot_size = 4,
                   out_height = 9, out_width = 7.5, alpha_values = rep(1, 3),
                   print_colormodel = print_colormodel)
coef_plot_function(start_ggplot = ggplot(res.s, aes(x=coef, y=v3, 
                                                    shape=v2, color=v2, alpha = v2)), 
                   file_name = paste0("Figure_23"),
                   main_title_name = NULL,
                   x_title_name = "Standardized Difference (Dem - NonDem)",
                   dot_size = 4,
                   out_height = 9, out_width = 7.5, alpha_values = rep(1, 3),
                   print_colormodel = print_colormodel)
# coefficient plot; ordinal spending variables; including military spending (standardized)
coef_plot_function(start_ggplot = ggplot(res.s.o, aes(x=coef, y=v3, 
                                                      shape=v2, color=v2, alpha = v2)), 
                   file_name = 
                     paste0("Figure_24"),
                   main_title_name = NULL, dot_size = 3.5,
                   x_title_name = "Standardized Difference (Dem - NonDem)",
                   out_height = 10, out_width = 7, 
                   alpha_values = rep(1, 3),
                   print_colormodel = print_colormodel)

# main cofficient plot (non-standardizd)
res.n <- res.s
res.n[,c(1,2)] <- data.frame(rbind(region.n, gdp.n, religion.n, white.n,
                                   oil.n, allies.n, trade.n, 
                                   exercise.n, invest.n, force.n))
coef_plot_function(start_ggplot = ggplot(res.n, aes(x=coef, y=v3, 
                                                    shape=v2, color=v2, alpha = v2)), 
                   file_name = 
                     paste0("Figure_25"),
                   main_title_name = NULL, scale_type = "free_x", 
                   x_title_name = "Non-standardized Difference (Dem - NonDem)",
                   dot_size = 3.5,
                   out_height = 10.5, out_width = 7, 
                   alpha_values = rep(1, 3),
                   print_colormodel = print_colormodel)

##################################
# make density plots
plot_data <- data.frame(Z = d$Z, V = d$V, d$region, d$gdp, d$religion, 
                        d$white, d$oil,
                        d$allies, d$trade, d$exercise, d$invest, d$force)
ncol(plot_data)-2

plot_data <- reshape2::melt(plot_data, measure.vars = 
                              names(plot_data)[3:ncol(plot_data)])
plot_data$variable <- relab(old.var = plot_data$variable, 
                            old.labs = unique(plot_data$variable),
                            new.labs = p.labs_l, text = TRUE)
all_combo <- expand.grid(Z = unique(plot_data$Z), V = unique(plot_data$V),
                         variable = unique(plot_data$variable))
# GDP
myden(plot_data = plot_data, 
      combo = all_combo[all_combo$variable == p.labs_l[2],], 
      bandwidth = 1,
      short_name = "Figure_13", xlab = "GDP per Capita (Log USD)",
      print_colormodel = print_colormodel)

# Religion
# relabel 
probint <- c(mean(0:20), mean(21:40), mean(41:60), mean(61:80), mean(81:100))
likelylab <- c("Very Unlikely", "Unlikely", "Odds About Even", 
               "Likely", "Very Likely")
# graphing distribution
myden(plot_data = plot_data,  
      combo = all_combo[all_combo$variable == p.labs_l[3],], 
      bandwidth = 20,
      short_name = "Figure_14", xlab = "Likelihood of Being Majority Christian",
      x.limits=c(0, 100), x.breaks=probint, x.labels=likelylab, 
      print_colormodel = print_colormodel)

# White
myden(plot_data = plot_data,  
      combo = all_combo[all_combo$variable == p.labs_l[4],], 
      bandwidth = 20,
      short_name = "Figure_16", xlab = "Likelihood of Being Majority White",
      x.limits=c(0, 100), x.breaks=probint, x.labels=likelylab,
      print_colormodel = print_colormodel)

# Oil
myden(plot_data = plot_data,  
      combo = all_combo[all_combo$variable == p.labs_l[5],], 
      bandwidth = 20,
      short_name = "Figure_15", 
      xlab = "Likelihood of Not Having Large Oil Reserves",
      x.limits=c(0, 100), x.breaks=probint, x.labels=likelylab,
      print_colormodel = print_colormodel)

# Military Allies
myden(plot_data = plot_data,  
      combo = all_combo[all_combo$variable == p.labs_l[6],], 
      bandwidth = 20,
      short_name = "Figure_17", 
      xlab = "Likelihood of Military Alliance with the U.S Since WW II",
      x.limits=c(0, 100), x.breaks=probint, x.labels=likelylab,
      print_colormodel = print_colormodel)

# Trade
myden(plot_data = plot_data,  
      combo = all_combo[all_combo$variable == p.labs_l[7],], 
      bandwidth = 2,
      short_name = "Figure_18", 
      xlab = "Total Trade Flow (Log Millions USD)",
      print_colormodel = print_colormodel)

# Joint Military Exercises
myden(plot_data = plot_data,  
      combo = all_combo[all_combo$variable == p.labs_l[8],], 
      bandwidth = 20,
      short_name = "Figure_19", 
      xlab = "Likelihood of Joint Military Exercise with U.S.",
      x.limits=c(0, 100), x.breaks=probint, x.labels=likelylab,
      print_colormodel = print_colormodel)

# FDI
sizelab <- c("Very Low", "Low", "Medium", "High", "Very High")
myden(plot_data = plot_data,  
      combo = all_combo[all_combo$variable == p.labs_l[9],], 
      bandwidth = 1,
      short_name = "Figure_20", 
      xlab = "Level of Investment in U.S. Businesses",
      x.limits=c(0, 6), x.breaks=1:5, x.labels=sizelab,
      print_colormodel = print_colormodel)

# Force
myden(plot_data = plot_data,  
      combo = all_combo[all_combo$variable == p.labs_l[10],], 
      bandwidth = 1.5,
      short_name = "Figure_21", xlab = "Spending on Military (Log Millions USD)",
      print_colormodel = print_colormodel)


######################################
# NPC
# Warning: This block of code takes a very long time to run.
# Uncomment to run.
# d$Z_num <- NA
# d$Z_num[d$Z == "Democracy"] <- 1
# d$Z_num[d$Z == "Non-democracy"] <- 0
# d <- as.data.frame(d)
# NPC_func <- function(vignette_type, comb_func, seed = 30294,
#                      test_stat = "DiffMeans", na_rm = TRUE) {
#   temp <- NPC(data = d[d$V == vignette_type,], tr.var = "Z_num", tr.label = 1,
#               n.perms = 10000,
#               comb.fun = comb_func, seed = seed, test.statistic = test_stat,
#               y.vars = c("region.s", "gdp.s", "religion.s", "oil.s",
#                          "white.s", "allies.s", "trade.s",
#                          "exercise.s", "invest.s"),
#               na.rm = na_rm)
#   return(temp[[1]])
# }
# 
# # check for missingness
# miniumcf_check <- lapply(X = levels(d$V), NPC_func, comb_func = "MinimumCF",
#                          test_stat = "DiffSumObs", na_rm = FALSE)
# normalcf_check <- lapply(X = levels(d$V), NPC_func, comb_func = "NormalCF",
#                        test_stat = "DiffSumObs", na_rm = FALSE)
# productcf_check <- lapply(X = levels(d$V), NPC_func, comb_func = "ProductCF",
#                         test_stat = "DiffSumObs", na_rm = FALSE)
# # NPC: DiffSumWithNA
# miniumcf_res1 <- lapply(X = levels(d$V), NPC_func, comb_func = "MinimumCF",
#                         test_stat = "DiffSumWithNA", na_rm = FALSE)
# normalcf_res1 <- lapply(X = levels(d$V), NPC_func, comb_func = "NormalCF",
#                         test_stat = "DiffSumWithNA", na_rm = FALSE)
# productcf_res1 <- lapply(X = levels(d$V), NPC_func, comb_func = "ProductCF",
#                          test_stat = "DiffSumWithNA", na_rm = FALSE)
# rbind(c(miniumcf_res2[[1]][19], miniumcf_res2[[2]][19], miniumcf_res2[[3]][19]),
#       c(normalcf_res2[[1]][19], normalcf_res2[[2]][19], normalcf_res2[[3]][19]),
#       c(productcf_res2[[1]][19], productcf_res2[[2]][19], productcf_res2[[3]][19]))
# # NPC: Difference of Means
# miniumcf_res2 <- lapply(X = levels(d$V), NPC_func, comb_func = "MinimumCF",
#                         test_stat = "DiffMeans", na_rm = TRUE)
# normalcf_res2 <- lapply(X = levels(d$V), NPC_func, comb_func = "NormalCF",
#                         test_stat = "DiffMeans", na_rm = TRUE)
# productcf_res2 <- lapply(X = levels(d$V), NPC_func, comb_func = "ProductCF",
#                          test_stat = "DiffMeans", na_rm = TRUE)
# npc2 <- rbind(
#       c(productcf_res2[[1]][19], productcf_res2[[2]][19], productcf_res2[[3]][19]),
#       c(normalcf_res2[[1]][19], normalcf_res2[[2]][19], normalcf_res2[[3]][19]),
#       c(miniumcf_res2[[1]][19], miniumcf_res2[[2]][19], miniumcf_res2[[3]][19]))

#######################################
# Treatment Measure 1 of Democracy
# full democracy
d$R1_1 <- multiq(v1 = d$dem_a_1, v2 = d$dem_b_1, v3 = d$dem_c_1)
d$R1_1 <- relab(old.var = d$R1_1, old.labs = 1:5, new.labs = probint)
# democracy
d$R1_2 <- multiq(d$dem_a_2, d$dem_b_2, d$dem_c_2)
d$R1_2 <- relab(old.var = d$R1_2, old.labs = 1:5, new.labs = probint)
# semi-democracy
d$R1_3 <- multiq(d$dem_a_3, d$dem_b_3, d$dem_c_3)
d$R1_3 <- relab(old.var = d$R1_3, old.labs = 1:5, new.labs = probint)
# semi-autocracy
d$R1_4 <- multiq(d$dem_a_4, d$dem_b_4, d$dem_c_4)
d$R1_4 <- relab(old.var = d$R1_4, old.labs = 1:5, new.labs = probint)
# full autocracy
d$R1_5 <- multiq(d$dem_a_5, d$dem_b_5, d$dem_c_5)
d$R1_5 <- relab(old.var = d$R1_5, old.labs = 1:5, new.labs = probint)
# normalize
for (i in 1:nrow(d)) {
  s <- d$R1_1[i] + d$R1_2[i] + d$R1_3[i] + d$R1_4[i] + d$R1_5[i]
  d$R1_1[i] <- d$R1_1[i]/s
  d$R1_2[i] <- d$R1_2[i]/s
  d$R1_3[i] <- d$R1_3[i]/s
  d$R1_4[i] <- d$R1_4[i]/s
  d$R1_5[i] <- d$R1_5[i]/s
}

regime_labs <- c("Fully Democratic", "Democratic", 
                 "Somewhat Democratic/\nSomewhat Non-democratic",
                 "Non-democratic", "Fully Non-democratic")
R1_summary <- d %>% group_by(V, Z) %>%
  dplyr::summarise(s_R1_1 = mean(R1_1, na.rm = TRUE),
            s_R1_2 = mean(R1_2, na.rm = TRUE),
            s_R1_3 = mean(R1_3, na.rm = TRUE),
            s_R1_4 = mean(R1_4, na.rm = TRUE),
            s_R1_5 = mean(R1_5, na.rm = TRUE)) %>% as.data.frame() %>%
  reshape2::melt()
R1_summary$variable <- relab(R1_summary$variable, 
                             unique(R1_summary$variable),
                             regime_labs, TRUE)
R1_summary$Z <- factor(as.character(R1_summary$Z), rev(levels(R1_summary$Z)))
R1_summary$variable <- factor(R1_summary$variable, 
                              levels = unique(R1_summary$variable))
# generate bar chart
ggplot(R1_summary, mapping = aes(y = value, x = variable, fill = Z)) + 
  geom_bar(stat = "identity", position = "dodge") + 
  facet_wrap( ~ V, ncol=2) + coord_flip() + theme_bb() + 
  ylab("Mean Probability") + xlab("Regime Type (Ordered by Polity Score)") +
  scale_fill_manual(values = c("steelblue4", "orange1"), name = "Treatment Assignment Regime Type") + ylim(0, 1) + 
  guides(fill = guide_legend(reverse = TRUE))
ggsave(paste0(images_directory, "Figure_26.pdf"), 
       width = 8, height = 6, dpi = 300, colormodel = print_colormodel)

# regression by regime type
regime_labs <- c("Fully Democratic", "Democratic", 
                 "Somewhat Democratic/Somewhat Non-democratic",
                 "Non-democratic", "Fully Non-democratic")
r_R1_1 <- myreg(Y = d$R1_1, Z = d$Z, V = d$V)
r_R1_2 <- myreg(Y = d$R1_2, Z = d$Z, V = d$V)
r_R1_3 <- myreg(Y = d$R1_3, Z = d$Z, V = d$V)
r_R1_4 <- myreg(Y = d$R1_4, Z = d$Z, V = d$V)
r_R1_5 <- myreg(Y = d$R1_5, Z = d$Z, V = d$V)
r_R1 <- data.frame(rep(regime_labs, each=3),
                   rep(levels(d$V), 5),
                   rbind(r_R1_1, r_R1_2, r_R1_3, r_R1_4, r_R1_5))
names(r_R1) <- c("placebo", "v", "coef", "se")
r_R1$placebo <- factor(as.character(r_R1$placebo), levels = rev(regime_labs))
r_R1$v2 <- factor(r_R1$v, levels = rev(levels(r_R1$v)))
r_R1_gg <- ggplot(r_R1, aes(x=coef, y=v2, 
                            shape=v, color=v, alpha = v))
coef_plot_function(start_ggplot = r_R1_gg, 
                   file_name = paste0("Figure_27"),
                   main_title_name = NULL,
                   x_title_name = "Difference in Likelihood (Dem - NonDem)",
                   dot_size = 3.5,
                   out_height = 6, out_width = 7, alpha_values = c(1,1,1),
                   print_colormodel = print_colormodel)

# check the probabilities sum up to 1
d$R1_sum <- d$R1_1 + d$R1_2 + d$R1_3 + d$R1_4 + d$R1_5
# impute Polity score
d$R1 <- 10*d$R1_1+mean(6:9)*d$R1_2+mean(1:5)*d$R1_3+
  mean(-5:0)*d$R1_4+mean(-10:-6)*d$R1_5
R1_sum <- d %>% group_by(V, Z) %>% dplyr::summarise(mean_R1 = mean(R1, na.rm = TRUE))
R1_sum$V <- factor(R1_sum$V, levels = rev(levels(R1_sum$V)))
# Make barchart
ggplot(R1_sum, aes(x = V, fill = Z, y = mean_R1)) +
  geom_bar(stat = "identity", position = "dodge") + coord_flip() +
  theme_bw() + theme(legend.position = "bottom") +
  scale_fill_manual(values = c("orange1", "steelblue4"), 
                    name = "Treatment Assignment Regime Type") +
  xlab("Vignette Type") + ylab("Imputed Polity Score")
ggsave(paste0(images_directory, "Figure_28_a.pdf"), 
       width = 8, height = 5, dpi = 300, colormodel = print_colormodel)
d$R1.s <- vscale(var = d$R1, vig = d$V)
# regression 
R1 <- myreg(Y = d$R1, Z = d$Z, V = d$V)
colnames(R1) <- c("coef", "SE")
R1_p <- data.frame(V = levels(d$V), R1)
names(R1_p)[2:3] <- c("coef", "se")
R1_p$V2 <- factor(R1_p$V, levels = rev(levels(R1_p$V)))
R1_coef <- ggplot(R1_p, mapping = aes(x = coef, y = V2, shape = V, color = V))
# make cofficient plot
coef_plot_function(start_ggplot = R1_coef, 
                   file_name = "Figure_28_b",
                   main_title_name = NULL, 
                   factet = FALSE,
                   x_title_name = "Difference in Imputed Polity Scores
                   (Dem - NonDem)",
                   out_height = 3, out_width = 5, alpha_values = c(1,1,1),
                   print_colormodel = print_colormodel)

# Treatment Measure 2: Characteristics of Democracy
R2_labs <- c("Elected Government",
             "Allows Opposition Party",
             "Free, Independent Media",
             "Freedom of Religion",
             "Has Executive Constraint",
             "Freedom of Assembly")
d$dc1 <- ifelse(d$demchar_1 == 1, TRUE, FALSE)
d$dc2 <- ifelse(d$demchar_4 == 1, TRUE, FALSE)
d$dc3 <- ifelse(d$demchar_5 == 1, TRUE, FALSE)
d$dc4 <- ifelse(d$demchar_6 == 1, TRUE, FALSE)
d$dc5 <- ifelse(d$demchar_7 == 1, TRUE, FALSE)
d$dc6 <- ifelse(d$demchar_8 == 1, TRUE, FALSE)
dc_sum <- d %>% group_by(Z, V) %>% 
  dplyr::summarize(dc1 = mean(dc1, na.rm = TRUE),
                   dc2 = mean(dc2, na.rm = TRUE),
                   dc3 = mean(dc3, na.rm = TRUE),
                   dc4 = mean(dc4, na.rm = TRUE),
                   dc5 = mean(dc5, na.rm = TRUE),
                   dc6 = mean(dc6, na.rm = TRUE)) %>% as.data.frame() %>%
  reshape2::melt()
dc_sum$variable <- relab(dc_sum$variable, old.labs = unique(dc_sum$variable),
                         new.labs = R2_labs, text = TRUE)
dc_sum$Z <- factor(dc_sum$Z, levels = rev(levels(dc_sum$Z)))
# run regressions
R2_res <- data.frame(v = rep(levels(d$V), 6),
                     placebo = rep(R2_labs, each = 3),
                     rbind(myreg(Y = d$dc1, Z = d$Z, V = d$V),
                           myreg(Y = d$dc2, Z = d$Z, V = d$V),
                           myreg(Y = d$dc3, Z = d$Z, V = d$V),
                           myreg(Y = d$dc4, Z = d$Z, V = d$V),
                           myreg(Y = d$dc5, Z = d$Z, V = d$V),
                           myreg(Y = d$dc6, Z = d$Z, V = d$V)))
names(R2_res)[3:4] <- c("coef", "se")
R2_res$v2 <- factor(R2_res$v, levels = rev(levels(R2_res$v)))
R2_gg <- ggplot(data = R2_res, mapping = aes(x = coef, y = v2, shape = v, color = v))
# make coefficient plot 
coef_plot_function(start_ggplot = R2_gg, 
                   file_name = "Figure_29_a",
                   main_title_name = NULL,
                   x_title_name = "Difference in Proportions Who Selected the Characteristic (Dem - NonDem)",
                   out_height = 6.5, out_width = 6, alpha_values = c(1,1,1),
                   print_colormodel = print_colormodel)

d$R2 <- d$dc1 + d$dc2 + d$dc3 + d$dc4 + d$dc5 + d$dc6
R2_count <- d %>% group_by(Z, V) %>% dplyr::summarize(count = mean(R2, na.rm = TRUE))
R2_count$V <- factor(R2_count$V, levels = rev(levels(R2_count$V)))
R2_count$Z <- factor(R2_count$Z, levels = rev(levels(R2_count$Z)))
d$R2.s <- vscale(var = d$R2, vig = d$V)
# regression
R2 <- myreg(Y = d$R2, Z = d$Z, V = d$V)
colnames(R2) <- c("coef", "SE")

R2_c <- data.frame(v = levels(d$V),
                   myreg(Y = d$R2, Z = d$Z, V = d$V))
names(R2_c)[2:3] <- c("coef", "se")
R2_c$v2 <- factor(R2_c$v, levels = rev(levels(R2_c$v)))
# make coefficient plot 
coef_plot_function(start_ggplot = ggplot(data = R2_c, 
                                         mapping = aes(x = coef, y = v2, 
                                                       shape = v, color = v)),
                   file_name = "Figure_29_b",
                   main_title_name = NULL, factet = FALSE,
                   x_title_name = "Difference in Number of Characteristics Selected (Dem - NonDem)",
                   out_height = 2, out_width = 5, alpha_values = c(1,1,1),
                   print_colormodel = print_colormodel)
# generate binary measure of democracy
d$bi_R <- ifelse(d$R1_1 + d$R1_2 > d$R1_4 + d$R1_5, TRUE, FALSE)
sum_bi_R <- d %>% group_by(Z, V) %>% dplyr::summarize(bi_R = mean(bi_R, na.rm = TRUE))
sum_bi_R$V <- factor(sum_bi_R$V, levels = rev(levels(sum_bi_R$V)))
sum_bi_R$Z <- factor(sum_bi_R$Z, levels = rev(levels(sum_bi_R$Z)))


# regression analysis
R1 <- myreg(Y = d$R1, Z = d$Z, V = d$V)
colnames(R1) <- c("coef", "SE")
bi_R <- myreg(Y = d$bi_R, Z = d$Z, V = d$V)
colnames(bi_R) <- c("coef", "SE")
combo_R <- rbind(bi_R, R1)
combo_R <- data.frame(V = levels(d$V), combo_R)
names(combo_R)[2:3] <- c("coef", "se")
combo_R$V2 <- factor(combo_R$V, levels = rev(levels(combo_R$V)))
combo_R$version <- c(rep("Dichotomous Regime Type", each = 3),
                     rep("Imputed Polity Score", each = 3))
combo_coef <- ggplot(combo_R, mapping = aes(x = coef, y = V2, shape = V, color = V))
f <- coef_plot_function(start_ggplot = combo_coef, print_out = FALSE,
                   main_title_name = NULL, factet_yes = FALSE,
                   x_title_name = "Difference in Proportions (Top) and Imputed Polity Scores (Bottom)\n(Dem-NonDem)", file_name = "combined_treatment_measure",
                   out_height = 4, out_width = 5, alpha_values = c(1,1,1),
                   print_colormodel = print_colormodel)
f + facet_wrap( ~ version, ncol=1, scales = "free_x") 
ggsave(paste0(images_directory, "Figure_08.pdf"),
       height = 3.5, width = 7, dpi = 300, colormodel = print_colormodel)

###############################
# ITT / IV Using R1
# ITT estimate of democracy on support for military action
d$support <- ifelse(is.na(d$support_1), d$support_2, d$support_1)
# coding #1: set "don't know" to NA
d$support1 <- d$support
d$support1[d$support>5] <- NA
d$support1 <- relab(old.var = d$support1, old.labs = 1:5, new.labs = 4:0)
support1 <- myreg(Y = d$support1, Z = d$Z, V = d$V)
colnames(support1) <- c("coef", "SE")
e_support1 <- myreg(Y = d$support1, Z = d$Z, V = d$ENE)
# coding #2: set "don't know" to 2
d$support2 <- d$support
d$support2[d$support>5] <- 2
d$support2 <- relab(old.var = d$support2, old.labs = 1:5, new.labs = 4:0)
support2 <- myreg(Y = d$support2, Z = d$Z, V = d$V)
colnames(support2)  <- c("coef", "SE")
e_support2 <- myreg(Y = d$support2, Z = d$Z, V = d$ENE)
# coding #3: binary coding for support
d$support3 <- d$support
d$support3 <- ifelse(d$support < 3, TRUE, FALSE)
support3 <- myreg(Y = d$support3, Z = d$Z, V = d$V)
e_support3 <- myreg(Y = d$support3, Z = d$Z, V = d$ENE)
colnames(support3)  <- c("coef", "SE")
if (output_print){
  print(support1)
  print(support2)
  print(support3)
}
###############################
# IV estimate of democracy on support for military action
iv_function <- function(support_DV, treatment_measure, V = d$V){
  IV_data <- na.omit(data.frame(support1 = support_DV, 
                                D = treatment_measure, Z = d$Z, V = V))
  iv.res1 <- myreg.iv(Y = IV_data$support1, D = IV_data$D,
                      Z = IV_data$Z, V = IV_data$V)
  colnames(iv.res1) <- c("coef", "SE")
  return(iv.res1)
}

iv.res1 <- iv_function(support_DV = d$support1*10, treatment_measure = d$R1)
iv.res2 <- iv_function(support_DV = d$support2*10, treatment_measure = d$R1)
iv.res3 <- iv_function(support_DV = d$support3*10, treatment_measure = d$R1)

iv.res1_b <- iv_function(support_DV = d$support1, treatment_measure = d$bi_R)
iv.res2_b <- iv_function(support_DV = d$support2, treatment_measure = d$bi_R)
iv.res3_b <- iv_function(support_DV = d$support3, treatment_measure = d$bi_R)

iv.res1_e <- iv_function(support_DV = d$support1*10, treatment_measure = d$R1,
                         V = d$ENE)
iv.res2_e <- iv_function(support_DV = d$support2*10, treatment_measure = d$R1,
                         V = d$ENE)
iv.res3_e <- iv_function(support_DV = d$support3*10, treatment_measure = d$R1,
                         V = d$ENE)

iv.res1_e_b <- iv_function(support_DV = d$support1, treatment_measure = d$bi_R,
                           V = d$ENE)
iv.res2_e_b <- iv_function(support_DV = d$support2, treatment_measure = d$bi_R,
                           V = d$ENE)
iv.res3_e_b <- iv_function(support_DV = d$support3, treatment_measure = d$bi_R,
                           V = d$ENE)


itt_iv_labs <- c("ITT", 
                 "IV: Imputed Polity Score Treatment Measure\n(Perceived Increase of 10 Polity Points)",
                 "IV: Dichotomous Treatment Measure\n(Perceived Non-democracy to Perceived Democracy)")
itt_iv_labs <- factor(itt_iv_labs, levels = itt_iv_labs)

# 2: ordinal scale 
support.res_2 <- data.frame(placebo = rep(itt_iv_labs, each = 3),
                            v = rep(levels(d$V), 3),
                            rbind(support2, iv.res2, iv.res2_b))
names(support.res_2)[3:4] <- c("coef", "se")
support.res_2$v2 <- factor(as.character(support.res_2$v), 
                           levels = rev(levels(support.res_2$v)))
support_gg <- ggplot(support.res_2, mapping = aes(x = coef, y = v2, color = v, shape = v))
coef_plot_function(start_ggplot = support_gg, 
                   file_name = paste0("Figure_30"),
                   main_title_name = 
                     "DV: Support for Using Force (Ordinal Scale 0 to 4)",
                   x_title_name = "Change in Support Using Force (Dem - NonDem)",
                   out_height = 4.5, out_width = 6, alpha_values = c(1,1,1),
                   print_colormodel = print_colormodel)

# 3: dichotomous measure
support.res_3 <- data.frame(placebo = rep(itt_iv_labs, each = 3),
                            v = rep(levels(d$V), 3),
                            rbind(support3, iv.res3, iv.res3_b))
names(support.res_3)[3:4] <- c("coef", "se")
support.res_3$v2 <- factor(as.character(support.res_3$v), 
                           levels = rev(levels(support.res_3$v)))
support_gg <- ggplot(support.res_3, mapping = 
                       aes(x = coef, y = v2, color = v, shape = v))
coef_plot_function(start_ggplot = support_gg, 
                   file_name = paste0("Figure_31"),
                   main_title_name = NULL,
                   x_title_name = 
                     "Change in Proportion Who Support Using Force (Dem - NonDem)",
                   out_height = 4.5, out_width = 6, alpha_values = c(1,1,1),
                   print_colormodel = print_colormodel)
coef_plot_function(start_ggplot = support_gg, 
                   file_name = paste0("Figure_03"),
                   main_title_name = NULL,
                   x_title_name = 
                     "Change in Proportion Who Support Using Force (Dem - NonDem)",
                   out_height = 4.5, out_width = 6, alpha_values = c(1,1,1),
                   print_colormodel = print_colormodel)
coef_plot_function(start_ggplot = support_gg, 
                   file_name = paste0("Figure_09"),
                   main_title_name = "DV: Support for Using Force (Dichotomous Measure)",
                   x_title_name = 
                     "Change in Proportion Who Support Using Force (Dem - NonDem)",
                   out_height = 4.5, out_width = 6, alpha_values = c(1,1,1),
                   print_colormodel = print_colormodel)

# different ENEs
# 2: ordinal scale
support.res_2_e <- data.frame(placebo = rep(itt_iv_labs, each = 2),
                              v = rep(levels(d$ENE), 3),
                              rbind(e_support2, iv.res2_e, iv.res2_e_b))
names(support.res_2_e)[3:4] <- c("coef", "se")
support.res_2_e$v2 <- factor(as.character(support.res_2_e$v), 
                             levels = rev(levels(support.res_2_e$v)))
support_gg_2_e <- ggplot(support.res_2_e, mapping = aes(x = coef, y = v2, color = v, shape = v))
coef_plot_function(start_ggplot = support_gg_2_e, 
                   file_name = paste0("Figure_32_a"),
                   main_title_name = 
                     "DV: Support for Using Force (Ordinal Scale 0 to 4)",
                   x_title_name = "Change in Support Using Force (Dem - NonDem)",
                   out_height = 4, out_width = 5.5, alpha_values = c(1,1,1),
                   print_colormodel = print_colormodel)
# 3: dichotomous measure
support.res_3_e <- data.frame(placebo = rep(itt_iv_labs, each = 2),
                              v = rep(levels(d$ENE), 3),
                              rbind(e_support3, iv.res3_e, iv.res3_e_b))
names(support.res_3_e)[3:4] <- c("coef", "se")
support.res_3_e$v2 <- factor(as.character(support.res_3_e$v), 
                             levels = rev(levels(support.res_3_e$v)))
support_gg_3_e <- ggplot(support.res_3_e, mapping = aes(x = coef, y = v2, color = v, shape = v))
coef_plot_function(start_ggplot = support_gg_3_e, 
                   file_name = paste0("Figure_32_b"),
                   main_title_name = 
                     "DV: Support for Using Force (Dichotomous Measure)",
                   x_title_name = "Change in Support Using Force (Dem - NonDem)",
                   out_height = 4, out_width = 5.5, alpha_values = c(1,1,1),
                   print_colormodel = print_colormodel)


############################################
# Balance Tests
# clean up the data
d$college <- ifelse(d$educ > 3, 1, 0)
d$democrat <- ifelse(d$partyid < 3, 1, 0)
d$age_num <- ifelse(d$age=="older than 100", 101, d$age+18)
d$male <- ifelse(d$sex == 1, 1, 0)
d$poliid_scale <- ifelse(d$poliid_1 > 0, d$poliid_1, NA)
# make the data for the bar chart
d_sum <- d %>% group_by(V, Z) %>%
  dplyr::summarise(age = mean(age_num, na.rm = TRUE),
            male = mean(male, na.rm = TRUE)*100,
            college = mean(college, na.rm = TRUE)*100,
            democrat = mean(democrat, na.rm = TRUE)*100,
            ideology = mean(poliid_scale, na.rm = TRUE)) %>% 
  as.data.frame() %>% reshape2::melt()
dem_var <- c("Mean Age", "Percent Male", "Percent with College Degree",
             "Percent Democrat", 
             "Mean Ideology Score\n(1 = Extremely Liberal; 7 = Extremely Conservative)")
d_sum$variable <- relab(d_sum$variable, old.labs = unique(d_sum$variable), 
                        new.labs = dem_var, text = TRUE)
d_sum$variable <- factor(d_sum$variable, levels = unique(d_sum$variable))
# make the bar chart
ggplot(data = d_sum, mapping = aes(x = V, y = value, fill = Z)) + 
  geom_bar(stat = "identity", position = "dodge") +
  facet_wrap(~ variable, ncol = 1, scale = "free_y") + theme_bb() + 
  xlab("Vignette Type") +  ylab("Value") +
  scale_fill_manual(values = c("orange1", "steelblue4"), name = "Treatment Assignment Regime Type")
ggsave(paste0(images_directory, "Figure_10.pdf"), 
       width = 6, height = 9, dpi = 300, colormodel = print_colormodel)

# balance tests
balance_res <- as.data.frame(matrix(NA, 15, 4))
names(balance_res) <- c("v", "placebo", "coef", "se")
balance_res$v <- rep(levels(d$V), each = 5)
balance_res$placebo <- factor(rep(dem_var, 3), levels = dem_var)
balance_F <- matrix(NA, 3, 3)
for (i in 1:3) {
  balance_res[1+(i-1)*5,3:4] <- robustse(lm(as.numeric(Z)-1 ~ age_num, 
                                            data = d[d$V==levels(d$V)[i],]))[2,1:2]
  balance_res[2+(i-1)*5,3:4] <- robustse(lm(as.numeric(Z)-1 ~ male, 
                                            data = d[d$V==levels(d$V)[i],]))[2,1:2]
  balance_res[3+(i-1)*5,3:4] <- robustse(lm(as.numeric(Z)-1 ~ college, 
                                            data = d[d$V==levels(d$V)[i],]))[2,1:2]
  balance_res[4+(i-1)*5,3:4] <- robustse(lm(as.numeric(Z)-1 ~ democrat, 
                                            data = d[d$V==levels(d$V)[i],]))[2,1:2]
  balance_res[5+(i-1)*5,3:4] <- robustse(lm(as.numeric(Z)-1 ~ poliid_scale, 
                                            data = d[d$V==levels(d$V)[i],]))[2,1:2] 
  overall <- felm(as.numeric(Z)-1 ~ college + democrat + age_num +
                    male + poliid_1, data = d[d$V==levels(d$V)[i],])
  b_overall <- summary(object = overall, robust = TRUE)$F
  balance_F[i,] <- c(levels(d$V)[i], 
                     paste0("F(", b_overall[2], ",", b_overall[3], ") = ", 
                            round(b_overall[1], 2)), round(b_overall[4], 3))
}
balance_res$v <- factor(balance_res$v, levels = unique(balance_res$v))
balance_res$v2 <- factor(as.character(balance_res$v), 
                         levels = rev(levels(balance_res$v)))
# make cofficient plot
balance_gg <- ggplot(balance_res, mapping = 
                       aes(x = coef, y = v2, color = v, shape = v))
coef_plot_function(start_ggplot = balance_gg, 
                   file_name = paste0("Figure_11"),
                   main_title_name = NULL,
                   x_title_name = 
                     "Difference of Means (Dem - NonDem)", scale_type = "free_x",
                   out_height = 7.5, out_width = 5.5, alpha_values = c(1,1,1),
                   print_colormodel = print_colormodel)
# make a regression table
# Table 2
stargazer(balance_F, summary = FALSE, rownames = FALSE)
                  
# further robustness checks
#####################################
# Abstract Encouragement Design
head(d)
d$abstract <- ifelse(!is.na(d$direct1), "Yes", NA)
d$abstract[!is.na(d$direct2)] <- "No"
d$abstract <- factor(d$abstract, levels = c("Yes", "No"))

std.abs <- data.frame(v = rep(levels(d$abstract), each = 10),
                      placebo = rep(p.labs_l, 2),
                      coef = NA, se = NA)
std.abs$v2 <- factor(std.abs$v, levels = rev(levels(std.abs$v)))

for (i in 1:2){
  std.abs[10*(i-1)+1, 3:4] <- robustse(lm(region.s ~ Z, 
                                          data = d[d$abstract == levels(d$abstract)[i],]))[2, 1:2]
  std.abs[10*(i-1)+2, 3:4] <- robustse(lm(gdp.s ~ Z, 
                                          data = d[d$abstract == levels(d$abstract)[i],]))[2, 1:2]
  std.abs[10*(i-1)+3, 3:4] <- robustse(lm(religion.s ~ Z, 
                                          data = d[d$abstract == levels(d$abstract)[i],]))[2, 1:2]
  std.abs[10*(i-1)+4, 3:4] <- robustse(lm(white.s ~ Z, 
                                          data = d[d$abstract == levels(d$abstract)[i],]))[2, 1:2]
  std.abs[10*(i-1)+5, 3:4] <- robustse(lm(oil.s ~ Z, 
                                          data = d[d$abstract == levels(d$abstract)[i],]))[2, 1:2]
  std.abs[10*(i-1)+6, 3:4] <- robustse(lm(allies.s ~ Z, 
                                          data = d[d$abstract == levels(d$abstract)[i],]))[2, 1:2]
  std.abs[10*(i-1)+7, 3:4] <- robustse(lm(trade.s ~ Z, 
                                          data = d[d$abstract == levels(d$abstract)[i],]))[2, 1:2]
  std.abs[10*(i-1)+8, 3:4] <- robustse(lm(exercise.s ~ Z, 
                                          data = d[d$abstract == levels(d$abstract)[i],]))[2, 1:2]
  std.abs[10*(i-1)+9, 3:4] <- robustse(lm(invest.s ~ Z,                                     
                                          data = d[d$abstract == levels(d$abstract)[i],]))[2, 1:2]
  std.abs[10*(i-1)+10, 3:4] <- robustse(lm(force.s ~ Z, 
                                           data = d[d$abstract == levels(d$abstract)[i],]))[2, 1:2]
}
coef_plot_function(start_ggplot = ggplot(std.abs, aes(x=coef, y=v, 
                                                      shape=v2, color=v2)), 
                   file_name = paste0("Supplementary_Appendix_Figure_3"),
                   main_title_name = NULL, 
                   color_values = c("purple", "orange"),
                   scale_name = "Abstract Encouragement Design",
                   x_title_name = "Standardized Difference  (Dem - NonDem)",
                   out_height = 8, out_width = 6, 
                   alpha_values = c(1,1,1),
                   print_colormodel = print_colormodel)
nonstd.abs <- std.abs
for (i in 1:2){
  nonstd.abs[10*(i-1)+1, 3:4] <- robustse(lm(region ~ Z, 
                                             data = d[d$abstract == levels(d$abstract)[i],]))[2, 1:2]
  nonstd.abs[10*(i-1)+2, 3:4] <- robustse(lm(gdp ~ Z, 
                                             data = d[d$abstract == levels(d$abstract)[i],]))[2, 1:2]
  nonstd.abs[10*(i-1)+3, 3:4] <- robustse(lm(religion ~ Z, 
                                             data = d[d$abstract == levels(d$abstract)[i],]))[2, 1:2]
  nonstd.abs[10*(i-1)+4, 3:4] <- robustse(lm(white ~ Z, 
                                             data = d[d$abstract == levels(d$abstract)[i],]))[2, 1:2]
  nonstd.abs[10*(i-1)+5, 3:4] <- robustse(lm(oil ~ Z, 
                                             data = d[d$abstract == levels(d$abstract)[i],]))[2, 1:2]
  nonstd.abs[10*(i-1)+6, 3:4] <- robustse(lm(allies ~ Z, 
                                             data = d[d$abstract == levels(d$abstract)[i],]))[2, 1:2]
  nonstd.abs[10*(i-1)+7, 3:4] <- robustse(lm(trade ~ Z, 
                                             data = d[d$abstract == levels(d$abstract)[i],]))[2, 1:2]
  nonstd.abs[10*(i-1)+8, 3:4] <- robustse(lm(exercise ~ Z, 
                                             data = d[d$abstract == levels(d$abstract)[i],]))[2, 1:2]
  nonstd.abs[10*(i-1)+9, 3:4] <- robustse(lm(invest ~ Z,                                     
                                             data = d[d$abstract == levels(d$abstract)[i],]))[2, 1:2]
  nonstd.abs[10*(i-1)+10, 3:4] <- robustse(lm(force ~ Z, 
                                              data = d[d$abstract == levels(d$abstract)[i],]))[2, 1:2]
}
# generate coefficient plot
coef_plot_function(start_ggplot = ggplot(nonstd.abs, aes(x=coef, y=v, 
                                                         shape=v2, color=v2)), 
                   file_name = paste0("Supplementary_Appendix_Figure_4"),
                   main_title_name = NULL, 
                   color_values = c("purple", "orange"),
                   scale_name = "Abstract Encouragement Design",
                   x_title_name = "Non-standardized Difference (Dem - NonDem)",
                   out_height = 9.5, out_width = 6, scale_type = "free_x",
                   alpha_values = c(1,1,1),
                   print_colormodel = print_colormodel)

# Berinsky Design: question block order
d$berinsky <- ifelse(d$DO_BL__A_Pre_outcomePlaceboTests != "", TRUE, FALSE)
unique(d$DO_BL__A_Pre_outcomePlaceboTests)
d$berinsky_order <- factor(d$DO_BL__A_Pre_outcomePlaceboTests)
berinsky_sum <- d %>% group_by(berinsky_order) %>% 
  dplyr::summarize(mean(support1, na.rm = TRUE)) %>% as.data.frame()
head(as.numeric(d$berinsky_order))
d$gdp_before <- rep(NA, nrow(d))
d$religion_before <- rep(NA, nrow(d))
d$regions_before <- rep(NA, nrow(d))
for (i in 1:nrow(d)){
  if (d$berinsky_order[i] != ""){
    d$gdp_before[i] <- gregexpr(pattern ='gdp', 
                                as.character(d$berinsky_order)[i])[[1]][1] <
      gregexpr(pattern ='support', as.character(d$berinsky_order)[i])[[1]][1]
    d$religion_before[i] <- gregexpr(pattern ='religion', 
                                     as.character(d$berinsky_order)[i])[[1]][1] <
      gregexpr(pattern ='support', as.character(d$berinsky_order)[i])[[1]][1]
    d$regions_before[i] <- gregexpr(pattern ='region', 
                                    as.character(d$berinsky_order)[i])[[1]][1] <
      gregexpr(pattern ='support', as.character(d$berinsky_order)[i])[[1]][1]
  }
}

d$before_num <- d$gdp_before + d$religion_before + d$regions_before
# Bivarate Regression
# support 1
b_labs <- factor(c("Regions", "GDP per Capita", "Religion"),
                 levels = c("Regions", "GDP per Capita", "Religion"))
b_sum1 <- data.frame(after = b_labs,
                     after = factor(b_labs, levels = rev(b_labs)),
                     rbind(robustse(lm(support1 ~ regions_before, data = d))[2,1:2],
                           robustse(lm(support1 ~ gdp_before, data = d))[2,1:2],
                           robustse(lm(support1 ~ religion_before, data = d))[2,1:2]))
names(b_sum1) <- c("after", "after1", "coef", "se")
# support 2
b_sum2 <- b_sum1
b_sum2[,3:4] <- rbind(robustse(lm(support2 ~ regions_before, data = d))[2,1:2],
                      robustse(lm(support2 ~ gdp_before, data = d))[2,1:2],
                      robustse(lm(support2 ~ religion_before, data = d))[2,1:2])
coef_plot_function(start_ggplot = ggplot(data = b_sum2, 
                                         mapping = aes(x = coef, y = after, shape = after, color = after)), 
                   file_name = "Supplementary_Appendix_Figure_6_a", 
                   x_title_name = "Difference in Support of Using Force",
                   factet_yes = FALSE, out_height = 2.5, out_width = 6, 
                   alpha_values = c(1,1,1), 
                   scale_name = "Placebo Test Question",
                   main_title_name = "DV: Support for Using Froce (Ordinal Scale: 0 to 4)",
                   print_colormodel = print_colormodel)
# support 3
b_sum3 <- b_sum2
b_sum3[,3:4] <- rbind(robustse(lm(support3 ~ regions_before, data = d))[2,1:2],
                      robustse(lm(support3 ~ gdp_before, data = d))[2,1:2],
                      robustse(lm(support3 ~ religion_before, data = d))[2,1:2])
coef_plot_function(start_ggplot = ggplot(data = b_sum3, 
                                         mapping = aes(x = coef, y = after, shape = after, color = after)), 
                   file_name = "Supplementary_Appendix_Figure_6_b", 
                   x_title_name = "Difference in Proportion Who Support of Using Force",
                   factet_yes = FALSE, out_height = 2.5, out_width = 6, 
                   alpha_values = c(1,1,1), 
                   scale_name = "Placebo Test Question",
                   main_title_name = "DV: Support for Using Force (Dichotomous Measure)",
                   print_colormodel = print_colormodel)

robustse(lm(support1 ~ gdp_before + religion_before + regions_before +
              gdp_before:religion_before + gdp_before:regions_before +
              religion_before:regions_before + 
              gdp_before:religion_before:regions_before, data = d))
robustse(lm(support1 ~ before_num, data = d))

all_2 <- robustse(lm(support2 ~ gdp_before + religion_before + 
                       regions_before +
                       gdp_before:religion_before + gdp_before:regions_before +
                       religion_before:regions_before + 
                       gdp_before:religion_before:regions_before, data = d)) 
robustse(lm(support2 ~ before_num, data = d))
robustse(lm(support3 ~ before_num, data = d))
all_3 <- robustse(lm(support3 ~ gdp_before + religion_before + regions_before +
                       gdp_before:religion_before + gdp_before:regions_before +
                       religion_before:regions_before + 
                       gdp_before:religion_before:regions_before, data = d)) 
stargazer(all_2, all_3, star.cutoffs = c(0.05, 0.01, 0.001))


# regular design
# placebo test block: FL_84
# outcome measure (support for war) block: FL_78
d$regular_placebo <- rep(NA, nrow(d))
for (i in 1:nrow(d)){
  if (d$DO_BR_FL_71[i] != "") {
    placebo <- gregexpr(pattern ='FL_84', as.character(d$DO_BR_FL_71)[i])[[1]][1] 
    support <- gregexpr(pattern ='FL_78', as.character(d$DO_BR_FL_71)[i])[[1]][1] 
    d$regular_placebo[i] = placebo < support 
  }
}
# regressions 
robustse(lm(support1 ~ regular_placebo, data = d))
robustse(lm(support2 ~ regular_placebo, data = d))
robustse(lm(support3 ~ regular_placebo, data = d))


interaction1_1 <- robustse(lm(support2 ~ Z + regions_before + Z:regions_before, data = d))
interaction1_2 <- robustse(lm(support2 ~ Z + gdp + Z:gdp, data = d))
interaction1_3 <- robustse(lm(support2 ~ Z + religion + Z:religion, data = d))
stargazer(interaction1_1, interaction1_2, interaction1_3, 
          dep.var.caption =  "Support for Using Force (Ordinal Scale)",
          star.cutoffs = c(0.05, 0.01, 0.001))

stargazer(
  robustse(lm(support3 ~ Z + regions_before + Z:regions_before, data = d)),
  robustse(lm(support3 ~ Z + gdp + Z:gdp, data = d)),
  robustse(lm(support3 ~ Z + religion + Z:religion, data = d)),
  dep.var.caption =  "DV: Support for Using Force (Dichotomous Measure)",
  star.cutoffs = c(0.05, 0.01, 0.001))

# Attention check
d$time_spent <- as.numeric((d$V9 - d$V8)/60)
outliers <- quantile(x = d$time_spent, probs = c(0.05, 0.95))
d$time_ok <- (d$time_spent > outliers[1] & d$time_spent < outliers[2])
d$a_passed <- ifelse(d$attention == 0, TRUE, FALSE)
d$c_order <- (ifelse(d$gdp_1 %in% c(1,7) | d$gdp_2 %in% c(1,7), TRUE, FALSE) +
                ifelse(d$religion_1 %in% c(1,5) | d$religion_2 %in% c(1,5), TRUE, FALSE) +
                ifelse(d$oil_1 %in% c(1,5) | d$oil_2 %in% c(1,5), TRUE, FALSE) +
                ifelse(d$white_1 %in% c(1,5) | d$white_2 %in% c(1,5), TRUE, FALSE) + 
                ifelse(d$alliance_1 %in% c(1,5) | d$alliance_2 %in% c(1,5), TRUE, FALSE) +
                ifelse(d$white_1 %in% c(1,5) | d$white_2 %in% c(1,5), TRUE, FALSE) +
                ifelse(d$trade_1 %in% c(1,5) | d$trade_2 %in% c(1,5), TRUE, FALSE) +
                ifelse(d$exercise_1 %in% c(1,5) | d$exercise_2 %in% c(1,5), TRUE, FALSE) +
                ifelse(d$invest_1 %in% c(1,5) | d$invest_2 %in% c(1,5), TRUE, FALSE) +
                ifelse(d$force_1 %in% c(1,5) | d$force_2 %in% c(1,5), TRUE, FALSE))/10
# Principal Axis Factor Analysis
# Pricipal Components Analysis
# entering raw data and extracting PCs
# from the correlation matrix
cor(data.frame(d$time_spent, d$a_passed, d$c_order), use = "complete.obs")
fit <- principal(data.frame(d$time_ok, d$a_passed, d$c_order),
                 nfactors=1, rotate="varimax")
fit # print results 
d$pca_score <- fit$score

robustse(lm(region.s ~ Z + pca_score + Z:pca_score, d))
robustse(lm(gdp.s ~ Z + pca_score + Z:pca_score, d))
robustse(lm(religion.s ~ Z + pca_score + Z:pca_score, d))
robustse(lm(white.s ~ Z + pca_score + Z:pca_score, d))
robustse(lm(oil.s ~ Z + pca_score + Z:pca_score, d))
robustse(lm(allies.s ~ Z + pca_score + Z:pca_score, d))
robustse(lm(trade.s ~ Z + pca_score + Z:pca_score, d))
robustse(lm(exercise.s ~ Z + pca_score + Z:pca_score, d))
robustse(lm(invest.s ~ Z + pca_score + Z:pca_score, d))
robustse(lm(force.s ~ Z + pca_score + Z:pca_score, d))

############################
# Question Order Effects

# Notes:
# Type A
# placebo test block: FL_84
# outcome measure (support for war) block: FL_78
# Type B
# placebo test block: FL_177
# outcome measure (support for war) block: FL_167

f_before <- function(variable, block_a, block_b) {
  res <- rep(NA, length(variable))
  for (i in 1:length(variable)) {
    res[i] <- gregexpr(pattern = block_a, as.character(variable)[i])[[1]][1] <
          gregexpr(pattern = block_b, as.character(variable)[i])[[1]][1]
  }
  return(res)
}

d$placebo_before <- ifelse(d$DO_BR_FL_173 == "", 
             f_before(d$DO_BR_FL_71, "FL_84", "FL_78"),
             NA)
d$placebo_after <- NA
d$placebo_after[d$placebo_before == TRUE] <- FALSE
d$placebo_after[d$placebo_before == FALSE] <- TRUE

d$plausibility_before <- ifelse(d$DO_BR_FL_173 == "", 
                                f_before(d$DO_BR_FL_71, "Plausibility check", "FL_78"),
                                NA)
d$manipulation_before <- ifelse(d$DO_BR_FL_173 == "", 
                                f_before(d$DO_BR_FL_71, "Manipulation Measure", "FL_78"),
                                NA)
d$open_write_before <- ifelse(d$DO_BR_FL_173 == "", 
                              f_before(d$DO_BR_FL_71, "Open Write", "FL_78"),
                              NA)
d$mediation_before <- ifelse(d$DO_BR_FL_78 == "", NA,
                             f_before(d$DO_BR_FL_78, "Mediation|Morality", "Support"))
# DV: Drop NA
order_1 <- rbind(robustse(lm(support1 ~ placebo_before, data = d))[2,1:2],
                 robustse(lm(support1 ~ open_write_before, data = d))[2,1:2],
  robustse(lm(support1 ~ manipulation_before, data = d)) [2,1:2],
  robustse(lm(support1 ~ plausibility_before, data = d))[2,1:2]) %>% 
  as.data.frame()
names(order_1) <- c("coef", "se")
block_labels <- c("Placebo Tests: Multiple Choice",
  "Placebo Test: Open-ended Response",
  "Treatment Measure",
  "Plausibility Check")
order_1$block <- factor(block_labels)
order_1$block2 <- factor(as.character(order_1$block),
                         levels = rev(levels(order_1$block)))
# DV: Missing as 2
order_2 <- rbind(robustse(lm(support2 ~ placebo_before, data = d[!d$berinsky,]))[2,1:2],
                 robustse(lm(support2 ~ open_write_before, data = d[!d$berinsky,]))[2,1:2],
                 robustse(lm(support2 ~ manipulation_before, data = d[!d$berinsky,])) [2,1:2],
                 robustse(lm(support2 ~ plausibility_before, data = d[!d$berinsky,]))[2,1:2]) %>% 
  as.data.frame()
names(order_2) <- c("coef", "se")
order_2$block <- factor(block_labels)
order_2$block2 <- factor(as.character(order_2$block),
                         levels = rev(levels(order_2$block)))
coef_plot_function(start_ggplot = ggplot(data = order_2, 
                                         mapping = aes(x = coef, y = block, shape = block2, color = block2)), 
                   file_name = "Supplementary_Appendix_Figure_5_a", 
                   stack_direction = "vertical",
                   legend_position = "top",
                   x_title_name = "Difference in Support of Using Force",
                   factet_yes = FALSE, out_height = 4, out_width = 6, 
                   alpha_values = c(1,1,1,1), 
                   color_values = 
                     c("firebrick3","forestgreen","dodgerblue3", "orange3"),
                   shape_values = c(21,22,23,24),
                   scale_name = "Block Before Support for Using Force Question",
                   main_title_name = "DV: Support for Using Force (Ordinal Scale: 0 to 4)",
                   print_colormodel = print_colormodel)
# DV: proportion
order_3 <- rbind(robustse(lm(support3 ~ placebo_before, data = d[!d$berinsky,]))[2,1:2],
                 robustse(lm(support3 ~ open_write_before, data = d[!d$berinsky,]))[2,1:2],
                 robustse(lm(support3 ~ manipulation_before, data = d[!d$berinsky,])) [2,1:2],
                 robustse(lm(support3 ~ plausibility_before, data = d[!d$berinsky,]))[2,1:2]) %>% 
  as.data.frame()
names(order_3) <- c("coef", "se")
order_3$block <- factor(block_labels)
order_3$block2 <- factor(as.character(order_3$block),
                         levels = rev(levels(order_3$block)))
coef_plot_function(start_ggplot = ggplot(data = order_3, 
                                         mapping = aes(x = coef, y = block, shape = block2, color = block2)), 
                   file_name = "Supplementary_Appendix_Figure_5_b", 
                   stack_direction = "vertical",
                   legend_position = "top",
                   x_title_name = "Difference in Proportion Who Support of Using Force",
                   factet_yes = FALSE, out_height = 4, out_width = 6, 
                   alpha_values = c(1,1,1,1), 
                   color_values = 
                     c("firebrick3","forestgreen","dodgerblue3", "orange3"),
                   shape_values = c(21,22,23,24),
                   scale_name = "Block Before Support for Using Force Question",
                   main_title_name = "DV: Support for Using Force (Dichotomous Measure)",
                   print_colormodel = print_colormodel)


robustse(lm(support1 ~ placebo_before + manipulation_before + open_write_before +
              placebo_before:manipulation_before + placebo_before:open_write_before +
              manipulation_before:open_write_before + 
              placebo_before:manipulation_before:open_write_before, data = d))

robustse(lm(support1 ~ mediation_before, data = d))
robustse(lm(support1 ~ Z + mediation_before + Z:mediation_before, data = d)) 

# check if support for war appears before placebo tests makes imbalance worse
placebo_vars <- c("region.s", "gdp.s", "religion.s", "oil.s",
  "white.s", "allies.s", "trade.s",
  "exercise.s", "invest.s", "force.s")
placebo_form <- paste0(placebo_vars, " ~ Z + placebo_after + Z:placebo_after")
order_imb <- data.frame(placebo = rep(p.labs_l, 3),
                        V = rep(levels(d$V), each = 10),
                        coef = NA,
                        se = NA)
for (i in 1:length(placebo_form)) {
  # Basic
  order_imb[i,3:4] <- robustse(lm(placebo_form[i], 
                               data = d[!d$berinsky & 
                                          d$V==levels(d$V)[1],]))[4,1:2]
  # Covariate Control
  order_imb[i+10,3:4] <- robustse(lm(placebo_form[i], 
                               data = d[!d$berinsky & 
                                          d$V==levels(d$V)[2],]))[4,1:2]
  # ENE
  order_imb[i+20,3:4] <- robustse(lm(placebo_form[i], 
                               data = d[!d$berinsky & 
                                          d$V==levels(d$V)[3],]))[4,1:2]
}
order_imb$V2 <- factor(order_imb$V, levels = rev(levels(order_imb$V)))
coef_plot_function(start_ggplot = ggplot(data = order_imb, 
                                         mapping = aes(x = coef, y = V2, shape = V, color = V)), 
                   file_name = "Supplementary_Appendix_Figure_7", 
                   legend_position = "top",
                   x_title_name = "Standardized Difference-in-Difference\n[(Dem - NonDem | Support for Force Question Before) - (Dem - NonDem | Support for Force Question After)]",
                   factet_yes = TRUE, out_height = 10, out_width = 9, 
                   scale_name = "Vignette Type",
                   main_title_name = NULL, alpha_values = c(1,1,1))

# print date and time of run
print("End date and time")
Sys.time()

